home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 7 / Apprentice-Release7.iso / Source Code / Pascal / Applications / NIH Image 1.62b11 / src / Background.p < prev    next >
Text File  |  1996-03-01  |  41KB  |  840 lines
  1. unit Background;
  2. {************************************************************************}
  3. {*     by Michael Castle and Janice Keller                                                                                                    *}
  4. {*     University of Michigan Mental Health Research Institute (MHRI)                                                     *}
  5. {*     e-mail address: mike.castle@umich.edu                                                                                  *}
  6. {************************************************************************}
  7. {*     Rolling ball and rolling disk algorithms inspired by Stanley Sternberg's article, "Biomedical        *}
  8. {*  Image Processing", IEEE Computer, January 1983.  Discussions with Rork Kuick and Tom               *}
  9. {*  Ford-Holevinski also contributed to the development of the algorithms and improvements in their   *}
  10. {*  efficiency.                                                                                                                                                *}
  11. {************************************************************************}
  12.  
  13. interface
  14.  
  15.     uses
  16.         Types, Memory, QuickDraw, Packages, Menus, Events, Fonts, Scrap, ToolUtils, Resources, Errors, Palettes, Windows,
  17.         globals, Utilities, Graphics, Camera, Filters;
  18.  
  19.     procedure DoBackgroundMenuEvent (MenuItem: integer);
  20.  
  21.  
  22. implementation
  23.  
  24.     type
  25.         IntRow = array[0..MaxLine] of Integer;
  26.         BackSubKindType = (RollingHorizontalArc, RollingVerticalArc, RollingBothArcs, RollingBall);
  27.  
  28.     var
  29.         ArcTrimPer: integer;                                  {trim off percentage of each side of the rolling ball patch}
  30.         shrinkfactor: integer;                                 {shrink the image and ball by this factor before rolling ball}
  31.         BackSubKind: BackSubKindType;                {which kind of background subtraction are we doing}
  32.         IntPlotWidth: Boolean;
  33.         Intplotwidthval: integer;
  34.         BoundRect: rect;
  35.         xxcenter, yycenter: integer;                                        {center of rectangular mask used in MinIn2DMask}
  36.         xmaskmin, ymaskmin: integer;                                   {upper left corner of mask used in AvgIn2DMask}
  37.         backgroundptr, ballptr, smallimageptr: ptr;              {ptrs to background, rolling ball, shrunk image memory}
  38.         backgroundaddr, balladdr, smallimageaddr: longint;   {addrs of background, rolling ball shrunk image memory}
  39.         patchwidth: LongInt;                                                     {x or y dimension of the rolling ball patch}
  40.         leftroll, rightroll, toproll, bottomroll: integer;         {bounds of the shrunk image}
  41.         Aborting: boolean;
  42.  
  43.  
  44.  
  45.  
  46.  
  47.     procedure RollBall;
  48. {*******************************************************************************}
  49. {*     RollBall 'rolls' a filtering object over a (shrunken) image in order to find the image's smooth continuous    *}
  50. {*  background.  For the purpose of explaining this algorithm, imagine that the 2D grayscale image has a third     *}
  51. {*  (height) dimension defined by the intensity value (0-255) at every point in the image.  The center of the      *}
  52. {*  filtering object, a patch from the top of a sphere having radius BallRadius, is moved along each scan line of     *}
  53. {*  the image so that the patch is tangent to the image at one or more points with every other point on the patch    *}
  54. {*  below the corresponding (x,y) point of the image.  Any point either on or below the patch during this process*}
  55. {*  is considered part of the background.  Shrinking the image before running this procedure is advised due to      *}
  56. {*  the fourth-degree complexity of the algorithm.  Care has been taken to avoid unnecessary operations (exp.      *}
  57. {*  multiplication inside loops) in this code.                                                                                                                *}
  58. {*******************************************************************************}
  59.         var
  60.             halfpatchwidth,                                {distance in x or y from patch center to any edge}
  61.             ptsbelowlastpatch,                           {number of points we may ignore because they were below last patch}
  62.             left, right, top, bottom,                   {}
  63.             xpt, ypt,                                           {current (x,y) point in the shrunken image}
  64.             xpt2, ypt2,                                      {current (x,y) point in the patch relative to upper left corner}
  65.             xval, yval,                                        {location in ball in shrunken image coordinates}
  66.             zdif,                                                  {difference in z (height) between point on ball and point on image}
  67.             zmin,                                                {smallest zdif for ball patch with center at current point}
  68.             zctr,                                                 {current height of the center of the sphere of which the patch is a part}
  69.             zadd: integer;                                    {height of a point on patch relative to the xy-plane of the shrunken image}
  70.             ballpt,                                               {index to chunk of memory storing the precomputed ball patch}
  71.             imgpt,                                               {index to chunk of memory storing the shrunken image}
  72.             backgrpt,                                          {index to chunk of memory storing the calculated background}
  73.             ybackgrpt,                                        {displacement to current background scan line}
  74.             ybackgrinc,                                      {distance in memory between two shrunken y-points in background}
  75.             smallimagewidth: longint;               {length of a scan line in shrunken image}
  76.             p1, p2: ptr;                                      {temporary pointers to background, ball, or small image}
  77.     begin
  78.         UpdateMeter(20, 'Finding Background...');
  79.         left := 1;
  80.         right := rightroll - leftroll - 1;
  81.         top := 1;
  82.         bottom := bottomroll - toproll - 1;
  83.         smallimagewidth := right - left + 3;
  84.         halfpatchwidth := patchwidth div 2;
  85.         ybackgrinc := shrinkfactor * (BoundRect.right - BoundRect.left);  {real dist btwn 2 adjacent (dy=1) shrunk pts}
  86.         zctr := 0;                                            {start z-center in the xy-plane}
  87.         for ypt := top to (bottom + patchwidth) do begin
  88.                 for xpt := left to (right + patchwidth) do {while patch is tangent to edges or within image...}
  89.                     begin                                           {xpt is far right edge of ball patch}
  90. {do we have to move the patch up or down to make it tangent to but not above image?...}
  91.                         zmin := 255;                              {highest could ever be 255}
  92.                         ballpt := balladdr;
  93.                         ypt2 := ypt - patchwidth;          {ypt2 is top edge of ball patch}
  94.                         imgpt := smallimageaddr + ypt2 * smallimagewidth + xpt - patchwidth;
  95.                         while ypt2 <= ypt do begin
  96.                                 xpt2 := xpt - patchwidth;      {xpt2 is far left edge of ball patch}
  97.                                 while xpt2 <= xpt do            {check every point on ball patch}
  98.                                     begin                                   {only examine points on }
  99.                                         if ((xpt2 >= left) and (xpt2 <= right) and (ypt2 >= top) and (ypt2 <= bottom)) then begin
  100.                                                 p1 := ptr(ballpt);
  101.                                                 p2 := ptr(imgpt);
  102.                                                 zdif := BAND(p2^, 255) - (zctr + BAND(p1^, 255));  {curve - circle points}
  103.                                                 if (zdif < zmin) then begin {keep most negative, since ball should always be below curve}
  104.                                                         zmin := zdif;
  105.                                                     end;
  106.                                             end;  {if xpt2,ypt2}
  107.                                         ballpt := ballpt + 1;          {step thru the ball patch memory}
  108.                                         xpt2 := xpt2 + 1;
  109.                                         imgpt := imgpt + 1;
  110.                                     end;  {while xpt2 }
  111.                                 ypt2 := ypt2 + 1;
  112.                                 imgpt := imgpt - patchwidth - 1 + smallimagewidth;
  113.                             end;  {while ypt2}
  114.                         if (zmin <> 0) then
  115.                             zctr := zctr + zmin;                {move ball up or down if we find a new minimum}
  116.                         if (zmin < 0) then
  117.                             ptsbelowlastpatch := halfpatchwidth    {ignore left half of ball patch when dz < 0}
  118.                         else
  119.                             ptsbelowlastpatch := 0;
  120. {now compare every point on ball with background,  and keep highest number}
  121.                         yval := ypt - patchwidth;
  122.                         ypt2 := 0;
  123.                         ballpt := balladdr;
  124.                         ybackgrpt := backgroundaddr + (yval - top + 1) * ybackgrinc;
  125.                         while ypt2 <= patchwidth do begin
  126.                                 xval := xpt - patchwidth + ptsbelowlastpatch;
  127.                                 xpt2 := ptsbelowlastpatch;
  128.                                 ballpt := ballpt + ptsbelowlastpatch;
  129.                                 backgrpt := ybackgrpt + (xval - left + 1) * shrinkfactor;
  130.                                 while xpt2 <= patchwidth do begin     {for all the points in the ball patch}
  131.                                         if ((xval >= left) and (xval <= right) and (yval >= top) and (yval <= bottom)) then begin
  132.                                                 p1 := ptr(ballpt);
  133.                                                 zadd := zctr + BAND(p1^, 255);
  134.                                                 p1 := ptr(backgrpt);
  135.                                                 if (zadd > BAND(p1^, 255)) then  {keep largest adjustment}
  136.                                                     p1^ := zadd;
  137.                                             end;
  138.                                         ballpt := ballpt + 1;
  139.                                         xval := xval + 1;
  140.                                         xpt2 := xpt2 + 1;
  141.                                         backgrpt := backgrpt + shrinkfactor;     {move to next point in x}
  142.                                     end;  {while xpt2}
  143.                                 yval := yval + 1;
  144.                                 ypt2 := ypt2 + 1;
  145.                                 ybackgrpt := ybackgrpt + ybackgrinc;       {move to next point in y}
  146.                             end;  {while ypt2}
  147.                     end;  {for xpt }
  148.                 if ((ypt mod 5) = 0) or not FasterBackgroundSubtraction then begin
  149.                         UpdateMeter(20 + (ord4(ypt - top) * 70) div (bottom + patchwidth - top), 'Finding Background...');
  150.                         if CommandPeriod then begin
  151.                                 beep;
  152.                                 Aborting := true;
  153.                                 Exit(RollBall);
  154.                             end;
  155.                     end;
  156.             end;  {for ypt}
  157.     end;
  158.  
  159.  
  160.     function MinIn2DMask {(xmaskmin,ymaskmin: integer)}
  161.         : integer;
  162. {*******************************************************************************}
  163. {*     MinInMask finds the minimum pixel value in a shrinkfactor X shrinkfactor mask.                                           *}
  164. {*******************************************************************************}
  165.         var
  166.             i, j,                                           {loop indices to step through mask}
  167.             thispixel,                                  {value at current pixel in mask}
  168.             min,                                          {temporary minimum value in mask}
  169.             nextrowoffset: integer;             {distance in memory from end of mask in this row to beginning in next}
  170.             paddr: longint;                           {address of current mask pixel}
  171.             p: ptr;                                        {pointer to current pixel in mask}
  172.     begin
  173.         with info^ do begin
  174.                 min := 255;
  175.                 nextrowoffset := bytesperrow - shrinkfactor;
  176.                 paddr := ord4(PicBaseAddr) + ymaskmin * bytesperrow + xmaskmin;
  177.                 for j := 1 to shrinkfactor do begin
  178.                         for i := 1 to shrinkfactor do begin
  179.                                 p := ptr(paddr);
  180.                                 thispixel := BAND(p^, 255);
  181.                                 if (thispixel < min) then
  182.                                     min := thispixel;
  183.                                 paddr := paddr + 1;
  184.                             end;     {for i}
  185.                         paddr := paddr + nextrowoffset;
  186.                     end;     {for j}
  187.                 MinIn2DMask := min;
  188.             end; {with}
  189.     end;
  190.  
  191.  
  192.     procedure GetRollingBall;
  193. {******************************************************************************}
  194. {*     This procedure computes the location of each point on the rolling ball patch relative to the center of the     *}
  195. {*  sphere containing it.  The patch is located in the top half of this sphere.  The vertical axis of the sphere         *}
  196. {*  passes through the center of the patch.  The projection of the patch in the xy-plane below is a square.           *}
  197. {******************************************************************************}
  198.         var
  199.             rsquare,                                                                         {rolling ball radius squared}
  200.             xtrim,                                                                            {# of pixels trimmed off each end of ball to make patch}
  201.             xval, yval,                                                                     {x,y-values on patch relative to center of rolling ball}
  202.             smallballradius, diam,                                                  {radius and diameter of rolling ball}
  203.             temp,                                                                             {value must be >=0 to take square root}
  204.             halfpatchwidth: integer;                                                {distance in x or y from center of patch to any edge}
  205.             i,                                                                                    {index into rolling ball patch memory}
  206.             ballsize: LongInt;                                                           {size of rolling ball memory}
  207.             p: ptr;                                                                            {pointer to current point on the ball patch}
  208.     begin
  209.         smallballradius := ballradius div shrinkfactor;           {operate on small-sized image with small-sized ball}
  210.         if smallballradius < 1 then
  211.             smallballradius := 1;
  212.         rsquare := smallballradius * smallballradius;
  213.         diam := smallballradius * 2;
  214.         xtrim := (ArcTrimPer * diam) div 100;                      {only use a patch of the rolling ball}
  215.         patchwidth := diam - xtrim - xtrim;
  216.         halfpatchwidth := smallballradius - xtrim;                   {this is half the patch width}
  217.         ballsize := (patchwidth + 1) * (patchwidth + 1);
  218.         ballptr := NewPtrClear(ballsize);
  219.         p := ballptr;
  220.         for i := 0 to ballsize - 1 do begin
  221.                 xval := i mod (patchwidth + 1) - halfpatchwidth;
  222.                 yval := i div (patchwidth + 1) - halfpatchwidth;
  223.                 temp := rsquare - (xval * xval) - (yval * yval);
  224.                 if (temp >= 0) then
  225.                     p^ := round(sqrt(temp))
  226.                 else
  227.                     p^ := 0;
  228.                 p := ptr(ord4(p) + 1);
  229.             end;
  230.     end;
  231.  
  232.  
  233.     procedure InterpolateBackground2D; {(leftroll, rightroll, toproll, bottomroll: integer; backgroundaddr: longint)}
  234. {******************************************************************************}
  235. {*     This procedure uses bilinear interpolation to find the points in the full-scale background given the points *}
  236. {*  from the shrunken image background.  Since the shrunken background is found from an image composed of    *}
  237. {*  minima (over a sufficiently large mask), it is certain that no point in the full-scale interpolated                 *}
  238. {*  background has a higher pixel value than the corresponding point in the original image.                                  *}
  239. {******************************************************************************}
  240.         var
  241.             i, ii,                                                   {horizontal loop indices}
  242.             j, jj,                                                  {vertical loop indices}
  243.             hloc, vloc,                                          {position of current pixel in calculated background}
  244.             vinc,                                                   {memory offset from current calculated pos to current interpolated pos}
  245.             lastvalue, nextvalue: integer;           {calculated pixel values between which we are interpolating}
  246.             p,                                                        {pointer to current interpolated pixel value}
  247.             bglastptr, bgnextptr: ptr;                 {pointers to calculated pixel values between which we are interpolating}
  248.             width: LongInt;
  249.     begin
  250.         vloc := 0;
  251.         with BoundRect do begin
  252.                 width := right - left;
  253.                 for j := 1 to bottomroll - toproll - 1 do begin     {interpolate to find background interior}
  254.                         hloc := 0;
  255.                         vloc := vloc + shrinkfactor;
  256.                         for i := 1 to rightroll - leftroll do begin
  257.                                 hloc := hloc + shrinkfactor;
  258.                                 bgnextptr := ptr(backgroundaddr + vloc * width + hloc);
  259.                                 bglastptr := ptr(ord4(bgnextptr) - shrinkfactor);
  260.                                 nextvalue := BAND(bgnextptr^, 255);
  261.                                 lastvalue := BAND(bglastptr^, 255);
  262.                                 for ii := 1 to shrinkfactor - 1 do begin     {interpolate horizontally}
  263.                                         p := ptr(ord4(bgnextptr) - ii);
  264.                                         p^ := lastvalue + (shrinkfactor - ii) * (nextvalue - lastvalue) div shrinkfactor;
  265.                                     end;     {for ii}
  266.                                 for ii := 0 to shrinkfactor - 1 do begin     {interpolate vertically}
  267.                                         bglastptr := ptr(backgroundaddr + (vloc - shrinkfactor) * width + hloc - ii);
  268.                                         bgnextptr := ptr(backgroundaddr + vloc * width + hloc - ii);
  269.                                         lastvalue := BAND(bglastptr^, 255);
  270.                                         nextvalue := BAND(bgnextptr^, 255);
  271.                                         vinc := 0;
  272.                                         for jj := 1 to shrinkfactor - 1 do begin
  273.                                                 vinc := vinc - right + left;
  274.                                                 p := ptr(ord4(bgnextptr) + vinc);
  275.                                                 p^ := lastvalue + (shrinkfactor - jj) * (nextvalue - lastvalue) div shrinkfactor;
  276.                                             end;     {for jj}
  277.                                     end;     {for ii}
  278.                             end;     {for i}
  279.                     end;     {for j}
  280.             end;   {with boundrect}
  281.     end;
  282.  
  283.  
  284.     procedure ExtrapolateBackground2D; {(leftroll, rightroll, toproll, bottomroll: integer; backgroundaddr: longint)}
  285. {******************************************************************************}
  286. {*     This procedure uses linear extrapolation to find pixel values on the top, left, right, and bottom edges of      *}
  287. {*  the background.  First it finds the top and bottom edge points by extrapolating from the edges of the                *}
  288. {*  calculated and interpolated background interior.  Then it uses the edge points on the new calculated,               *}
  289. {*  interpolated, and extrapolated background to find all of the left and right edge points.  If extrapolation yields *}
  290. {*  values below zero or above 255, then they are set to zero and 255 respectively.                                              *}
  291. {******************************************************************************}
  292.         var
  293.             ii, jj,                                                 {horizontal and vertical loop indices}
  294.             hloc, vloc,                                          {position of current pixel in calculated/interpolated background}
  295.             edgeslope,                                          {difference of last two consecutive pixel values on an edge}
  296.             pvalue,                                               {current extrapolated pixel value}
  297.             lastvalue, nextvalue: integer;           {calculated pixel values from which we are extrapolating}
  298.             p,                                                        {pointer to current extrapolated pixel value}
  299.             bglastptr, bgnextptr: ptr;                 {pointers to calculated pixel values from which we are extrapolating}
  300.             width: LongInt;
  301.     begin
  302.         with BoundRect do begin
  303.                 width := right - left;
  304.                 for hloc := shrinkfactor to shrinkfactor * (rightroll - leftroll) - 1 do begin     {extrapolate on top and bottom}
  305.                         bglastptr := ptr(backgroundaddr + shrinkfactor * width + hloc);
  306.                         bgnextptr := ptr(backgroundaddr + (shrinkfactor + 1) * width + hloc);
  307.                         lastvalue := BAND(bglastptr^, 255);
  308.                         nextvalue := BAND(bgnextptr^, 255);
  309.                         edgeslope := nextvalue - lastvalue;
  310.                         p := bglastptr;
  311.                         pvalue := lastvalue;
  312.                         for jj := 1 to shrinkfactor do begin
  313.                                 p := ptr(ord4(p) - right + left);
  314.                                 pvalue := pvalue - edgeslope;
  315.                                 if (pvalue < 0) then
  316.                                     p^ := 0
  317.                                 else if (pvalue > 255) then
  318.                                     p^ := 255
  319.                                 else
  320.                                     p^ := pvalue;
  321.                             end;     {for jj}
  322.                         bglastptr := ptr(backgroundaddr + (shrinkfactor * (bottomroll - toproll - 1) - 1) * width + hloc);
  323.                         bgnextptr := ptr(backgroundaddr + shrinkfactor * (bottomroll - toproll - 1) * width + hloc);
  324.                         lastvalue := BAND(bglastptr^, 255);
  325.                         nextvalue := BAND(bgnextptr^, 255);
  326.                         edgeslope := nextvalue - lastvalue;
  327.                         p := bgnextptr;
  328.                         pvalue := nextvalue;
  329.                         for jj := 1 to (bottom - top - 1) - shrinkfactor * (bottomroll - toproll - 1) do begin
  330.                                 p := ptr(ord4(p) + right - left);
  331.                                 pvalue := pvalue + edgeslope;
  332.                                 if (pvalue < 0) then
  333.                                     p^ := 0
  334.                                 else if (pvalue > 255) then
  335.                                     p^ := 255
  336.                                 else
  337.                                     p^ := pvalue;
  338.                             end;     {for jj}
  339.                     end;     {for hloc}
  340.                 for vloc := top to bottom - 1 do begin     {extrapolate on left and right}
  341.                         bglastptr := ptr(backgroundaddr + (vloc - top) * width + shrinkfactor);
  342.                         bgnextptr := ptr(ord4(bglastptr) + 1);
  343.                         lastvalue := BAND(bglastptr^, 255);
  344.                         nextvalue := BAND(bgnextptr^, 255);
  345.                         edgeslope := nextvalue - lastvalue;
  346.                         p := bglastptr;
  347.                         pvalue := lastvalue;
  348.                         for ii := 1 to shrinkfactor do begin
  349.                                 p := ptr(ord4(p) - 1);
  350.                                 pvalue := pvalue - edgeslope;
  351.                                 if (pvalue < 0) then
  352.                                     p^ := 0
  353.                                 else if (pvalue > 255) then
  354.                                     p^ := 255
  355.                                 else
  356.                                     p^ := pvalue;
  357.                             end;     {for ii}
  358.                         bgnextptr := ptr(backgroundaddr + (vloc - top) * width + shrinkfactor * (rightroll - leftroll - 1) - 1);
  359.                         bglastptr := ptr(ord4(bgnextptr) - 1);
  360.                         lastvalue := BAND(bglastptr^, 255);
  361.                         nextvalue := BAND(bgnextptr^, 255);
  362.                         edgeslope := nextvalue - lastvalue;
  363.                         p := bgnextptr;
  364.                         pvalue := nextvalue;
  365.                         for ii := 1 to (right - left - 1) - shrinkfactor * (rightroll - leftroll - 1) + 1 do begin
  366.                                 p := ptr(ord4(p) + 1);
  367.                                 pvalue := pvalue + edgeslope;
  368.                                 if (pvalue < 0) then
  369.                                     p^ := 0
  370.                                 else if (pvalue > 255) then
  371.                                     p^ := 255
  372.                                 else
  373.                                     p^ := pvalue;
  374.                             end;     {for ii}
  375.                     end;     {for vloc}
  376.             end;   {with BoundRect}
  377.     end;
  378.  
  379.  
  380.     procedure SubtractBackground2D;
  381. {*****************************************************************************}
  382. {*     This procedure subtracts each pixel from the calculated/interpolated/extrapolated background from the  *}
  383. {*  corresponding pixel value in the original image.  The resulting image is stored in place of the original        *}
  384. {*  image.  Any pixel subtractions with results below zero are given the value zero.                                           *}
  385. {*****************************************************************************}
  386.         var
  387.             hloc, vloc,                                          {current pixel location in image and background}
  388.             pvalue: integer;                                 {difference at current pixel location}
  389.             offset,                                                 {offset in memory from beginning of original image to current scan line}
  390.             backgrpt: LongInt;                              {offset to current point in background}
  391.             p: ptr;                                                {temporary pointer to image or background points}
  392.             Databand: Linetype;                           {current scan line in image}
  393.             ControlKey: boolean;
  394.     begin
  395.         backgrpt := 0;
  396.         ControlKey := ControlKeyDown;
  397.         with Info^, BoundRect do begin
  398.                 for vloc := top to bottom - 1 do begin
  399.                         GetLine(0, vloc, pixelsperline, Databand);
  400.                         for hloc := left to right - 1 do begin
  401.                                 p := ptr(backgroundaddr + backgrpt);
  402.                                 pvalue := Databand[hloc] - BAND(p^, 255);
  403.                                 if ControlKey then
  404.                                     pvalue := BAND(p^, 255);
  405.                                 if pvalue < 0 then
  406.                                     Databand[hloc] := 0
  407.                                 else
  408.                                     Databand[hloc] := pvalue;
  409.                                 backgrpt := backgrpt + 1;
  410.                             end;     {for}
  411.                         offset := vloc * BytesPerRow;
  412.                         p := ptr(ord4(PicBaseAddr) + offset);
  413.                         BlockMove(@Databand, p, pixelsperline);
  414.                     end;  {for}
  415.             end;     {with}
  416.     end;
  417.  
  418.  
  419.     procedure Background2D;
  420. {******************************************************************************}
  421. {*     This procedure implements a rolling-ball algorithm for the removal of smooth continuous background       *}
  422. {*  from a two-dimensional gel image.  It rolls the ball (actually a square patch on the top of a sphere) on a       *}
  423. {*  low-resolution (by a factor of 'shrinkfactor' times) copy of the original image in order to increase speed     *}
  424. {*  with little loss in accuracy.  It uses interpolation and extrapolation to blow the shrunk image to full size.     *}
  425. {******************************************************************************}
  426.         var
  427.             tport: Grafptr;
  428.             i,                                     {loop index for shrunk image memory}
  429.             backgroundsize,              {size of the background memory}
  430.             smallimagesize: LongInt;     {size of the shrunk image memory}
  431.             p: ptr;                             {pointer to current pixel in shrunk image memory}
  432.             table: FateTable;             {not used}
  433.             width: LongInt;
  434.     begin
  435.         ShowWatch;
  436.         UpdateMeter(0, 'Building Rolling Ball...');
  437.         GetPort(tPort);
  438.         with Info^ do begin
  439.                 SetPort(GrafPtr(osPort));
  440.                 BoundRect := roiRect;
  441.             end;
  442.         GetRollingBall;                                                                  {precompute the rolling ball}
  443.         UpdateMeter(3, 'Building Rolling Ball...');
  444.         balladdr := ord4(ballptr);
  445.         with BoundRect do begin
  446.                 width := right - left;
  447.                 leftroll := left div shrinkfactor;                                  {left and right edges of shrunken image or roi}
  448.                 rightroll := right div shrinkfactor - 1;                      {on which to roll ball}
  449.                 toproll := top div shrinkfactor;
  450.                 bottomroll := bottom div shrinkfactor - 1;
  451.                 backgroundsize := width * (bottom - top);
  452.                 backgroundptr := NewPtrClear(backgroundsize);
  453.                 Aborting := backgroundptr = nil;
  454.                 backgroundaddr := ord4(backgroundptr);
  455.                 smallimagesize := ord4(rightroll - leftroll + 1);
  456.                 smallimagesize := smallimagesize * (bottomroll - toproll + 1);
  457.                 smallimageptr := NewPtrClear(smallimagesize);
  458.                 Aborting := Aborting or (smallimageptr = nil);
  459.                 smallimageaddr := ord4(smallimageptr);
  460.                 if not aborting then begin
  461.                         UpdateMeter(6, 'Smoothing Image ');
  462.                         filter(unweightedAvg, 1, table);                                {smooth image before shrinking}
  463.                         UpdateMeter(10, concat('Shrinking Image ', long2str(shrinkfactor), ' times...'));
  464.                         for i := 0 to smallimagesize - 1 do begin                {create a lower resolution image for ball-rolling}
  465.                                 p := ptr(smallimageaddr + i);
  466.                                 xmaskmin := left + shrinkfactor * (i mod (rightroll - leftroll + 1));
  467.                                 ymaskmin := top + shrinkfactor * (i div (rightroll - leftroll + 1));
  468.                                 p^ := MinIn2DMask;                            {each point in small image is minimum of its neighborhood}
  469.                             end;
  470.                         if not aborting then begin
  471.                                 Undo;        {restore original unsmoothed image}
  472.                                 RollBall;
  473.                             end;
  474.                     end
  475.                 else
  476.                     beep;
  477.                 if not Aborting then begin
  478.                         UpdateMeter(90, 'Interpolating Background...');
  479.                         InterpolateBackground2D;                              {interpolate to find background interior}
  480.                         UpdateMeter(95, 'Extrapolating Background...');
  481.                         ExtrapolateBackground2D;                             {extrapolate on top and bottom}
  482.                         UpdateMeter(98, 'Subtracting Background...');
  483.                         SubtractBackground2D;                                  {subtract background from original image}
  484.                         UpdateMeter(100, 'Subtracting Background...');
  485.                     end;
  486.             end;   {with boundrect}
  487.         DisposePtr(backgroundptr);                           {free up background, rolling ball, shrunk image memory}
  488.         DisposePtr(ballptr);
  489.         DisposePtr(smallimageptr);
  490.         DisposeWindow(MeterWindow);
  491.         MeterWindow := nil;
  492.         SetPort(tPort);
  493.     end;
  494.  
  495.  
  496.     procedure RollArc (left, rightminusone, diam: integer; var background, ballpoints: IntRow; var Dataline: Linetype);
  497.         var
  498.             xpt, xpt2, xval, ydif, ymin, yctr, bpt, yadd: integer;
  499.     begin
  500.         for xpt := left to rightminusone do begin
  501.                 background[xpt] := -255;         {init background curve to minimum values}
  502.             end;
  503.         yctr := 0;                                   {start y-center at the x axis}
  504.         for xpt := left to (rightminusone + diam - 1) do {while semicircle is tangent to edges or within curve...}
  505.             begin                                       {xpt is far right edge of semi-circle}
  506. {do we have to move the circle?...}
  507.                 ymin := 256;                          {highest could ever be 255}
  508.                 bpt := 0;
  509.                 xpt2 := xpt - diam;                {xpt2 is far left edge of semi-circle}
  510.                 while xpt2 <= xpt do            {check every point on semicircle}
  511.                     begin
  512.                         if ((xpt2 >= left) and (xpt2 <= rightminusone)) then begin  {only examine points on curve}
  513.                                 ydif := dataline[xpt2] - (yctr + ballpoints[bpt]);                {curve minus circle points}
  514.                                 if (ydif < ymin) then begin {keep most negative, since ball should always be below curve}
  515.                                         ymin := ydif;
  516.                                     end;
  517.                             end;  {if xpt2 }
  518.                         bpt := bpt + 1;
  519.                         xpt2 := xpt2 + 1;
  520.                     end;  {while xpt2 }
  521.                 if (ymin <> 256) then{if we found a new minimum...}
  522.                     yctr := yctr + ymin;   {move circle up or down}
  523. {now compare every point on semi with background,  and keep highest number}
  524.                 xval := xpt - diam;
  525.                 xpt2 := 0;
  526.                 while xpt2 <= diam do begin  {for all the points in the semicircle}
  527.                         if ((xval >= left) and (xval <= rightminusone)) then begin
  528.                                 yadd := yctr + ballpoints[xpt2];
  529.                                 if (yadd > background[xval]) then  {keep largest adjustment}
  530.                                     background[xval] := yadd;
  531.                             end;
  532.                         xval := xval + 1;
  533.                         xpt2 := xpt2 + 1;
  534.                     end;  {while xpt2}
  535.             end;  {for xpt }
  536.     end;
  537.  
  538.  
  539.     function MinIn1DMask (var Databand: LineType; xcenter: integer): integer;
  540. {*******************************************************************************}
  541. {*     MinIn1DMask finds the minimum pixel value in a (2*shrinkfactor-1) mask about the point xcenter in the *}
  542. {*  current line.  This code must run FAST because it gets called OFTEN!                                                                   *}
  543. {*******************************************************************************}
  544.         var
  545.             i,                                                                              {index to pixels in the mask}
  546.             temp: integer;                                                          {temporary minimum value}
  547.     begin
  548.         temp := Databand[xcenter - shrinkfactor + 1];
  549.         for i := xcenter - shrinkfactor + 2 to xcenter + shrinkfactor - 1 do
  550.             if (Databand[i] < temp) then
  551.                 temp := Databand[i];
  552.         MinIn1DMask := temp;
  553.     end;
  554.  
  555.  
  556. {******************************************************************************}
  557. {*     This procedure computes the location of each point on the rolling arc relative to the center of the circle     *}
  558. {*  containing it.  The arc is located in the top half of this circle.  The vertical axis of the circle passes through  *}
  559. {*  the midpoint of the arc.  The projection of the arc on the x-axis below is a line segment.                                 *}
  560. {******************************************************************************}
  561.     procedure GetRollingArc (var arcpoints: IntRow; var arcwidth: integer);
  562.         var
  563.             xpt,                                                                                 {x-point along arc}
  564.             xval,                                                                               {x-point in arc array}
  565.             rsquare,                                                                         {shrunken arc radius squared}
  566.             xtrim,                                                                            {points to be trimmed off each end of arc}
  567.             smallballradius,                                                            {radius of shrunken arc which actually rolls}
  568.             diam: integer;                                                                 {diameter of shrunken arc's circle}
  569.     begin
  570.         smallballradius := ballradius div shrinkfactor;            { operate on small-sized image with small-sized ball}
  571.         rsquare := smallballradius * smallballradius;
  572.         for xpt := -smallballradius to smallballradius do        { find the ballpoints for arc based at  (x,y)=(0,0) }
  573.             begin
  574.                 xval := xpt + smallballradius;                                     {offset, can't have negative index}
  575.                 arcpoints[xval] := round(sqrt(abs(rsquare - (xpt * xpt))));  {Ys are positive, top half of circle}
  576.             end;
  577.         diam := smallballradius * 2;
  578.         xtrim := (ArcTrimPer * diam) div 100;                       {how many points to trim off each end}
  579.         arcwidth := diam - xtrim - xtrim;
  580.         for xpt := -smallballradius to smallballradius - xtrim - xtrim do begin
  581.                 xval := xpt + smallballradius;
  582.                 arcpoints[xval] := arcpoints[xval + xtrim];
  583.             end;
  584.         for xpt := smallballradius - xtrim - xtrim + 1 to smallballradius do begin
  585.                 xval := xpt + smallballradius;
  586.                 arcpoints[xval] := 0;
  587.             end;
  588.     end;
  589.  
  590.  
  591.     procedure ExtrapolateBackground1D (var Backline, Dataline: LineType; background: IntRow; leftroll, rightroll: integer);
  592. {******************************************************************************}
  593. {*     This procedure uses linear extrapolation to find pixel values on the left and right edges of the current        *}
  594. {*  line of the background.  It finds the edge points by extrapolating from the edges of the calculated and               *}
  595. {*  interpolated background interior.  If extrapolation yields values below zero or above 255, then they are set *}
  596. {*  to zero and 255 respectively.                                                                                                                               *}
  597. {******************************************************************************}
  598.         var
  599.             i,                                                                             {index to edges of background array}
  600.             hloc,                                                                       {}
  601.             edgeslope: integer;                                                 {}
  602.     begin
  603.         with BoundRect do begin
  604.                 edgeslope := (background[leftroll + 1] - background[leftroll + 2]) div shrinkfactor;
  605.                 for i := left to shrinkfactor * (leftroll + 1) - 1 do begin     {extrapolate on left edge}
  606.                         hloc := shrinkfactor * (leftroll + 1) - 1 + left - i;
  607.                         if (Backline[hloc + 1] + edgeslope < 0) then
  608.                             Backline[hloc] := 0
  609.                         else if (Backline[hloc + 1] + edgeslope > Dataline[hloc]) then
  610.                             Backline[hloc] := Dataline[hloc]
  611.                         else
  612.                             Backline[hloc] := Backline[hloc + 1] + edgeslope;
  613.                     end;
  614.                 edgeslope := (background[rightroll] - background[rightroll - 1]) div shrinkfactor;
  615.                 for hloc := shrinkfactor * rightroll to right - 1 do begin     {extrapolate on right edge}
  616.                         if (Backline[hloc - 1] + edgeslope < 0) then
  617.                             Backline[hloc] := 0
  618.                         else if (Backline[hloc - 1] + edgeslope > Dataline[hloc]) then
  619.                             Backline[hloc] := Dataline[hloc]
  620.                         else
  621.                             Backline[hloc] := Backline[hloc - 1] + edgeslope;
  622.                     end;
  623.             end;     {with}
  624.     end;
  625.  
  626.     procedure Background1D;
  627.         var
  628.             tport: Grafptr;
  629.             hloc, arcwidth, leftroll, rightroll, numpixels: integer;
  630.             left, right, top, bottom: integer;                      {image bounds; ROTATED if RollingVerticalArc}
  631.             i, j, maskwidth: integer;
  632.             background, arcpoints: IntRow;
  633.             vloc, offset: LongInt;
  634.             p: ptr;
  635.             Dataline, Backline, Smalldataline: Linetype;
  636.             str: str255;
  637.     begin
  638.         ShowWatch;
  639.         UpdateMeter(0, concat('Shrinking Image ', long2str(shrinkfactor), ' times...'));
  640.         GetPort(tPort);
  641.         with Info^ do begin
  642.                 SetPort(GrafPtr(osPort));
  643.                 BoundRect := roiRect;
  644.             end;
  645.         GetRollingArc(arcpoints, arcwidth);
  646.         maskwidth := shrinkfactor + shrinkfactor - 1;
  647.         case BackSubKind of
  648.             RollingHorizontalArc:  begin
  649.                     left := BoundRect.left;
  650.                     top := BoundRect.top;
  651.                     right := BoundRect.right;
  652.                     bottom := BoundRect.bottom;
  653.                     numpixels := Info^.pixelsperline;
  654.                     str := 'Rolling Disk Horizontally...';
  655.                 end;
  656.             RollingVerticalArc:  begin
  657.                     left := BoundRect.top;
  658.                     top := BoundRect.left;
  659.                     right := BoundRect.bottom;
  660.                     bottom := BoundRect.right;
  661.                     numpixels := Info^.nlines;
  662.                     str := 'Rolling Disk Vertically...';
  663.                 end;
  664.         end;     {case}
  665.         leftroll := left div shrinkfactor;                                  {left and right edges of shrunken image or roi}
  666.         rightroll := right div shrinkfactor - 1;                      {on which to roll arc}
  667.         for vloc := top to bottom - 1 do begin  {for ROI}
  668.                 case BackSubKind of
  669.                     RollingHorizontalArc: 
  670.                         GetLine(0, vloc, numpixels, Dataline);
  671.                     RollingVerticalArc: 
  672.                         GetColumn(vloc, 0, numpixels, Dataline);
  673.                 end;     {case}
  674.                 for i := leftroll + 1 to rightroll do
  675.                     smalldataline[i] := MinIn1DMask(Dataline, shrinkfactor * i - 1);
  676.                 RollArc(leftroll + 1, rightroll, arcwidth, background, arcpoints, smalldataline);  {roll arc on one line}
  677.                 for i := leftroll + 1 to rightroll do begin           {interpolate to find interior background points}
  678.                         hloc := shrinkfactor * i - 1;
  679.                         Backline[hloc] := background[i];
  680.                         for j := 1 to shrinkfactor - 1 do
  681.                             Backline[hloc - j] := background[i - 1] + (shrinkfactor - j) * (background[i] - background[i - 1]) div shrinkfactor;
  682.                     end;
  683.                 ExtrapolateBackground1D(Backline, Dataline, background, leftroll, rightroll);
  684.                 for i := left to right - 1 do begin                                {subtract background from current scan line}
  685.                         Dataline[i] := Dataline[i] - Backline[i];
  686.                         if Dataline[i] < 0 then
  687.                             Dataline[i] := 0;
  688.                     end;
  689.                 case BackSubKind of
  690.                     RollingHorizontalArc: 
  691.                         with Info^ do begin
  692.                                 offset := vloc * BytesPerRow;
  693.                                 p := ptr(ord4(PicBaseAddr) + offset);
  694.                                 BlockMove(@Dataline, p, numpixels);            {fast whole line write}
  695.                             end;
  696.                     RollingVerticalArc: 
  697.                         PutColumn(vloc, 0, numpixels, Dataline);         {slow whole column write}
  698.                 end;     {case}
  699.                 if ((vloc mod 8) = 0) and (vloc > 16) then begin
  700.                         UpdateMeter(((vloc - top) * 100) div (bottom - top - 1), str);
  701.                         if CommandPeriod then begin
  702.                                 beep;
  703.                                 Aborting := true;
  704.                                 leave;
  705.                             end;
  706.                     end;
  707.             end;
  708.         UpdateMeter(100, str);
  709.         DisposeWindow(MeterWindow);
  710.         MeterWindow := nil;
  711.         SetPort(tPort);
  712.     end;
  713.  
  714.     procedure SetUpGel;
  715.         var
  716.             frame: rect;
  717.             AutoSelectAll: boolean;
  718.             p: ptr;
  719.     begin
  720.         if NotinBounds or NotRectangular then
  721.             exit(SetUpGel);
  722.         StopDigitizing;
  723.         AutoSelectAll := not Info^.RoiShowing;
  724.         if AutoSelectAll then
  725.             SelectAll(false);
  726.         SetupUndoFromClip;
  727.         with info^ do begin
  728.                 frame := roiRect;
  729.                 if ((LutMode = GrayScale) or (LutMode = CustomGrayscale)) and (not IdentityFunction) then
  730.                     ApplyLookupTable;
  731.                 changes := true;
  732.             end;
  733.         case BackSubKind of
  734.             RollingHorizontalArc, RollingVerticalArc: 
  735.                 Background1D;    {--------------> call background subtract <-------------------}
  736.             RollingBall: 
  737.                 Background2D;
  738.             RollingBothArcs:  begin
  739.                     BackSubKind := RollingHorizontalArc;           {remove horizontal streaks}
  740.                     Background1D;
  741.                     BackSubKind := RollingVerticalArc;               {remove vertical streaks}
  742.                     if not aborting then
  743.                         Background1D;
  744.                     BackSubKind := RollingBothArcs;                   {leave BackSubKind as we found it}
  745.                 end;
  746.         end;     {case}
  747.         UpdatePicWindow;
  748.         SetUpRoiRect;
  749.         WhatToUndo := UndoFilter;
  750.         Info^.changes := true;
  751.         ShowWatch;
  752.         if AutoSelectAll then
  753.             KillRoi;
  754.         if Aborting then begin
  755.                 Undo;
  756.                 WhatToUndo := NothingToUndo;
  757.                 UpdatePicWindow;
  758.             end;
  759.     end;
  760.  
  761.  
  762.     procedure GetBallRadius;
  763.         var
  764.             SaveRadius: integer;
  765.             canceled: boolean;
  766.     begin
  767.         SaveRadius := BallRadius;
  768.         BallRadius := GetInt('Rolling BallRadius:', BallRadius, canceled);
  769.         if (BallRadius < 1) or (BallRadius > 319) or canceled then begin
  770.                 BallRadius := SaveRadius;
  771.                 if not canceled then
  772.                     beep;
  773.             end;
  774.     end;
  775.  
  776.  
  777.     procedure DoBackgroundMenuEvent (MenuItem: integer);
  778.         var
  779.             map_array: Ptr;
  780.     begin
  781.         if MenuItem <= RemoveStreaksItem then
  782.             if not CheckCalibration
  783.                 then exit(DoBackgroundMenuEvent);
  784.         MeterWindow := nil;
  785.         Aborting := false;
  786.         case MenuItem of
  787.             HorizontalItem, VerticalItem:  begin
  788.                     if FasterBackgroundSubtraction then begin
  789.                             ArcTrimPer := 20;
  790.                             shrinkfactor := 4;
  791.                         end
  792.                     else begin
  793.                             ArcTrimPer := 10;
  794.                             shrinkfactor := 2;
  795.                         end;
  796.                     if MenuItem = HorizontalItem then
  797.                         BackSubKind := RollingHorizontalArc
  798.                     else
  799.                         BackSubKind := RollingVerticalArc;
  800.                     SetUpGel;
  801.                 end;
  802.             Sub2DItem:  begin
  803.                     if FasterBackgroundSubtraction then begin
  804.                             if BallRadius > 15 then begin
  805.                                     ArcTrimPer := 20;     {trim 40% in x and y}
  806.                                     shrinkfactor := 8;
  807.                                 end
  808.                             else begin
  809.                                     ArcTrimPer := 16;  {trim 32% in x and y}
  810.                                     shrinkfactor := 4;
  811.                                 end
  812.                         end
  813.                     else begin  {faster not checked}
  814.                             if BallRadius > 15 then begin
  815.                                     ArcTrimPer := 16;  {trim 32% in x and y}
  816.                                     shrinkfactor := 4;
  817.                                 end
  818.                             else begin
  819.                                     ArcTrimPer := 12;   {trim 24% in x and y}
  820.                                     ShrinkFactor := 2;
  821.                                 end
  822.                         end;
  823.                     BackSubKind := RollingBall;
  824.                     SetUpGel;
  825.                 end;
  826.             RemoveStreaksItem:  begin
  827.                     ArcTrimPer := 20;
  828.                     shrinkfactor := 4;
  829.                     BackSubKind := RollingBothArcs;
  830.                     SetUpGel;
  831.                 end;
  832.             FasterItem: 
  833.                 FasterBackgroundSubtraction := not FasterBackgroundSubtraction;
  834.             RadiusItem: 
  835.                 GetBallRadius;
  836.         end; {case}
  837.     end;
  838.  
  839.  
  840. end.